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ABSTRACT 

We examine aspects of primordial star formation in the presence of a molec- 
ular hydrogen-dissociating ultraviolet background. We compare a set of AMR 
hydrodynamic cosmological simulations using a single cosmological realization 
but with a range of ultraviolet background strengths in the Lyman- Werner band. 
This allows us to study the effects of Lyman- Werner radiation on suppressing 
H2 cooling at low densities as well as the high-density evolution of the collapsing 
cloud core in a self- consistent cosmological framework. We find that the addi- 
tion of a photodissociating background results in a delay of the collapse of high 
density gas at the center of the most massive halo in the simulation and, as a 
result, an increase in the virial mass of this halo at the onset of baryon collapse. 
We find that, contrary to previous results. Population III star formation is not 
suppressed for J21 > 0.1, but occurs even with backgrounds as high as J21 = 1. 
We find that H2 cooling leads to collapse despite the depressed core molecular 
hydrogen fractions due to the elevated H2 cooling rates at T = 2 — 5 x 10'^ K. We 
observe a relationship between the strength of the photodissociating background 
and the rate of accretion onto the evolving protostellar cloud core, with higher 
LW background fluxes resulting in higher accretion rates. Finally, we find that 
the collapsing cloud cores in our simulations do not fragment at densities be- 
low n 10^*^ cm~^ regardless of the strength of the LW background, suggesting 
that Population III stars forming in halos with T^r ~ 10'^ K may still form in 
isolation. 
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1. Introduction 

The study of the formation of Population III stars in a cosmological context via high- 
resolution simulations is becoming a mature discipline, with multiple groups finding similar 
results (Bromm, Coppi, & Larson 2002; Bromm & Loeb 2004; Abel, Bryan, & Norman 
2002; O'Shea & Norman 2007; Yoshida et al. 2006b; Gao et al. 2006). However, all of 
these calculations operate under a common fundamental assumption - namely, the absence 
of an ultraviolet background. Given the abundance of halos in which Population III stars 
would form and the general consensus that these stars must be massive, one expects that, 
though only a small fraction of the volume of the universe would be ionized, a significant 
background of ultraviolet radiation in the Lyman- Werner (LW) band (11.18-13.6 eV), which 
is capable of photodissociating molecular hydrogen, would be present for the formation of 
the bulk of Population III stars (Omukai 2001; Yoshida et al. 2003; Machacek, Bryan, & 
Abel 2001; Wise & Abel 2005). Given that atomic hydrogen is optically thin to this radiation 
and that it easily destroys molecular hydrogen, this could have a significant impact on the 
formation of primordial stars. Machacek et al. (2001) used cosmological AMR simulations 
to study the formation of primordial stars in the presence of a soft UV background, and 
later including an x-ray background (Machacek, Bryan, & Abel 2003). They found that 
the LW background delays the formation of Population III stars and shifts halo formation 
to higher masses. However, their calculations were not of sufficiently high resolution to 
give any significant information about the formation of the protostellar cloud at the core 
of each collapsing cosmological halo. Yoshida et al. (2003) performed a suite of large SPH 
simulations of high redshift structure formation, including two calculations with a soft UV 
background with J21 = 0.01 and 0.1, where J21 is the mean intensity of the UV background 
in the LW band in units of 10~^^ erg s~^ cm~^ Hz~^ sr~^. They found that gas coohng is 
suppressed in the higher LW background case, and predicted that star formation would not 
occur. On the other hand. Oh & Haiman (2002) showed analytically that primordial gas 
in halos with T^r > 10"^ K can still collapse to high densities in the presence of a strong 
molecular hydrogen-dissociating background via atomic hydrogen line radiation, and can 
eventually form H2. They referred to these objects as "second generation objects" in the 
sense that although still of primordial composition, their formation pathway is different from 
the first stars. 

An issue that needs to be addressed concerns the mode of primordial star formation 
in halos with T^j^ ~ lO'' K and above. It is apparent from recent studies using large 
numbers of simulations that smaller halos, with masses of ~ 10^ — 10^ M©, appear to 
always form primordial stars in isolation (Gao et al. 2006; O'Shea & Norman 2007). Is 
this still the case in halos whose masses are an order of magnitude or more larger? Does a 
single primordial star form, or several? This is particularly relevant given that many semi- 
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analytic models of structure formation that follow the chemical evolution of structures at 
high redshift (Scannapieco et al. 2005; Schneider et al. 2006) suggest that halos of this scale 
which have not been enriched by metals (and thus will form primordial stars) exist up to at 
least z ~ 5. If this is true, there is expected to be a strong photodissociating background at 
that time (Yoshida et al. 2003; Wise & Abel 2005). 

This paper is the second in a series. In O'Shea & Norman (2007) (hereafter referred to as 
Paper I), we examined several aspects of Population III star formation in a ACDM universe, 
in the absence of an ultraviolet background which photodissociates molecular hydrogen. 
In this paper, we study Population III star formation in a single cosmological realization 
varying the strength of the LW background, with the goal of investigating the effect that 
this background has on the evolution and properties of gas in the halo core (at radii <C 1 pc) 
and the ultimate fate of the gas at the center of T^j^ ~ 10^ K halos. We use a single 
cosmological realization, varying the strength of the photodissociating background over the 
range suggested by Wise & Abel (2005), and look for the cutoff in Population III star 
formation suggested by Yoshida et al. (2003). 

In agreement with the findings of (Machacek et al. 2001), we find that the addition 
of a photodissociating background results in a delay of the collapse of high density gas at 
the center of the most massive halo in the simulation and, as a result, an increase in the 
virial mass of this halo at the onset of baryon collapse. We find that, contrary to the results 
of Yoshida et al. (2003), star formation is not suppressed for J 21 > 0.1, but occurs even 
with backgrounds as high as J21 = 1. We find that H2 cooling leads to collapse despite 
the depressed halo core molecular hydrogen fractions ~ 10~^ by two multiplicative 
effects: (1) the elevated H2 cooling rates per molecule at T = 2 — 5 x 10^ K, and (2) time. 
We find that halo core collapse occurs in the usual way once the gas in the halo core has 
become dense enough that the cooling time becomes much less than the Hubble time. We 
also observe a relationship between the strength of the photodissociating background and 
the rate of accretion onto the evolving protostellar cloud core, with higher LW background 
fluxes resulting in higher accretion rates. This is a simple consequence of the suppression 
of molecular hydrogen formation (and thus suppression of cooling) by the photodissociating 
background, as well as the higher virial temperatures of the more massive halos at the 
epoch of collapse. This may have implications for the range of Population III stellar masses. 
Finally, we find that the collapsing halo cores in our simulations do not fragment at densities 
below n ~ 10^° cm~^ regardless of the strength of the soft UV background, suggesting that 
Population III stars forming in halos with T^r ~ 10*^ K may still form in isolation. 

The organization of this paper is as follows. In Section [2] we provide a description 
of Enzo, the code used to perform the calculations in this paper and of the simulation 
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setup. The results from our simulations are presented in Sections [3] through [3 Section [3] 
discusses some of the mean halo properties observed in the calculations, Section H] discusses 
the evolution of the halo core prior to collapse, Section O compares spherically-averaged halo 
properties for all simulations at the epoch of collapse. Section discusses a variety of halo 
properties at a fixed redshift, and Section [7] compares the evolution of two representative 
simulations taken from our ensemble. In Section [8] we discuss neglected physics and possible 
numerical issues, and in Section [9] we discuss some of the results presented in this work and 
their implications. Finally, we present a summary of the main results in Section (TUl 



2. Methodology 

2.1. The Enzo code 

'EnzoQ is a publicly available, extensively tested adaptive mesh refinement (AMR) cos- 
mology code developed by Greg Bryan and others (Bryan & Norman 1997a, 1997b; Norman 
& Bryan 1999; O'Shea et al. 2004, 2005b). The specifics of the Enzo code are described in 
detail in these papers (and references therein), but we present a brief description here for 
clarity. 

The Enzo code couples an N-body particle-mesh (PM) solver (Efstathiou et al. 1985; 
Hockney & Eastwood 1988) used to follow the evolution of a collisionless dark matter compo- 
nent with an Eulerian AMR method for ideal gas dynamics by Berger & Colella (1989), which 
allows high dynamic range in gravitational physics and hydrodynamics in an expanding uni- 
verse. This AMR method (referred to as structured AMR) utilizes an adaptive hierarchy of 
grid patches at varying levels of resolution. Each rectangular grid patch (referred to as a 
"grid") covers some region of space in its parent grid which requires higher resolution, and 
can itself become the parent grid to an even more highly resolved child grid. Enzo's im- 
plementation of structured AMR places no fundamental restrictions on the number of grids 
at a given level of refinement, or on the number of levels of refinement. However, owing to 
limited computational resources it is practical to institute a maximum level of refinement, 
imax- Additionally, the Enzo AMR implementation allows arbitrary integer ratios of parent 
and child grid resolution, though in general for cosmological simulations (including the work 
described in this paper) a refinement ratio of 2 is used. 

Since the addition of more highly refined grids is adaptive, the conditions for refinement 
must be specified. In Enzo, the criteria for refinement can be set by the user to be a combi- 
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nation of any or all of the following: baryon or dark matter overdensity threshold, minimum 
resolution of the local Jeans length, local density gradients, local pressure gradients, local 
energy gradients, shocks, and cooling time. A cell reaching any or all of the user-specified 
criteria will then be flagged for refinement. Once all cells of a given level have been examined, 
rectangular solid boundaries are determined which minimally encompass the flagged cells on 
that level. A refined grid patch is then introduced within each such bounding volume, and 
the results are interpolated to a higher level of resolution. 

In Enzo, resolution of the equations being solved is adaptive in time as well as in space. 
The timcstcp in Enzo is satisfied on a Icvcl-by-lcvcl basis by finding the largest timcstcp 
such that the Courant condition (and an analogous condition for the dark matter particles) 
is satisfied by every cell on that level. All cells on a given level are advanced using the 
same timestep. Once a level L has been advanced in time Ai^, all grids at level L + 1 are 
advanced, using the same criteria for timestep calculations described above, until they reach 
the same physical time as the grids at level L. At this point grids at level L + 1 exchange 
baryon fiux information with their parent grids, providing a more accurate solution on level 
L. Cells at level L + 1 are then examined to see if they should be refined or de-refined, and 
the entire grid hierarchy is rebuilt at that level (including all more highly refined levels). 
The timestepping and hierarchy rebuilding processes are repeated recursively on every level 
to the maximum existing grid level in the simulation. 

Two different hydrodynamic methods are implemented in Enzo: the piecewisc parabolic 
method (PPM) (Woodward & Colclla 1984), which was extended to cosmology by Bryan 
et al. (1995), and the hydrodynamic method used in the ZEUS magnetohydrodynamics code 
(Stone & Norman 1992a, 1992b). We direct the interested reader to the papers describing 
both of these methods for more information, and note that PPM is the preferred choice 
of hydro method since it is higher-order-accurate and is based on a technique that does 
not require artificial viscosity, which smoothes shocks and can smear out features in the 
hydrodynamic flow. 

The chemical and cooling properties of primordial (metal-free) gas are followed using 
the method of Abel et al. (1997) and Anninos et al. (1997). This method follows the non- 
equilibrium evolution of a gas of primordial composition with 9 total species: H, H~^, He, 
i7e+, He^^ , H~, H2, H2, and e~. The code also calculates radiative heating and cooling 
following atomic line excitation, recombination, collisional excitation, free-free transitions, 
molecular line cooling, and Compton scattering of the cosmic microwave background, as well 
as any of approximately a dozen different models for a metagalactic ultraviolet background 
that heat the gas via photoionization and/or photodissociation. We model the cooling pro- 
cesses detailed in Abel et al. (1997), but use the GaUi & Palla (1998) molecular hydrogen 
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cooling function. The multispecies rate equations are solved out of equilibrium to properly 
model situations where, e.g., the cooling time of the gas is much shorter than the hydrogen 
recombination time. A total of 9 kinetic equations are solved, including 29 kinetic and radia- 
tive processes, for the 9 species mentioned above. The chemical reaction equation network 
is technically challenging to solve due to the huge range of reaction time scales involved; the 
characteristic creation and destruction time scales of the various species and reactions can 
differ by many orders of magnitude. As a result, the set of rate equations is extremely stiff, 
and an explicit scheme for integration of the rate equations can be costly if small enough 
timestep are taken to keep the network stable. This makes an implicit scheme preferable 
for such a set of equations, and Enzo solves the rate equations using a method based on a 
backwards differencing formula (BDF) in order to provide a stable and accurate solution. 

It is important to note the regime in which this chemistry model is valid. According 
to Abel et al. (1997) and Anninos et al. (1997), the reaction network is vahd for temper- 
atures between 10° — 10^ K. The original model discussed in these two references is only 
valid up to uh ~ 10^ cm"^. However, addition of the 3-body H2 formation process al- 
lows correct modeling of the gas chemistry up until the point where collisionally induced 
emission from molecular hydrogen becomes an important cooling processes, which occurs at 
uh ~ 10^"^ cm~^. We do not include heating by molecular hydrogen formation, which will be 
significant at densities of ~ 10^ cm~^ and above, and may affect temperature evolution at 
these high densities. A further concern is that the optically thin approximation for radiative 
coohng breaks down beginning at n^^ ~ 10^° — 10^^ cm~^. Beyond this point, modifications 
to the coohng function that take into account the non-negligible opacity of the gas to fine 
radiation from molecular hydrogen must be made, as discussed by Ripamonti & Abel (2004). 
Even with these modifications, a more correct description of the cooling of gas of primor- 
dial composition at high densities will require some form of radiation transport, which will 
greatly increase the cost of the simulations. 



2.2. Simulation setup 

The simulations discussed in this paper are set up in a similar way to those in Paper I. 
A dark matter-only calculation with 128^ particles in a three-dimensional simulation volume 
which is 0.6 h~^ Mpc (comoving) on a side is set up at z = 99 assuming a "concordance" 
cosmological model with no baryons: Qm = ^dm = 0.3, Qi, = 0.0, Q\ = 0.7, h = 0.7 (in 
units of 100 km/s/Mpc), as = 0.9, and using an Eisenstein & Hu power spectrum (Eiscnstein 
& Hu 1999) with a spectral index oi n — 1. The cold dark matter cosmological model is 
assumed. This calculation is then evolved to 2; = 15 using a maximum of four levels of 
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adaptive mesh refinement, refining on a dark matter overdensity of 8.0. At z = 15, tlie Hop 
fialo finding algoritlim (Eisenstein & Hut 1998) is used to find the most massive halo in the 
simulation. 

At this point, we generate a new set of initial conditions which contain the same large- 
scale power as the dark matter-only calculation, but include both dark matter and baryons 
such that the Lagrangian volume in which the halo in the dark matter-only calculation 
formed is resolved at high spatial and mass resolution using a series of static nested grids, 
with a 128^ root grid and three static nested grids, for an overall effective root grid size 
of 1024^ cells. The highest resolution grid is 256^ grid cells, and corresponds to a volume 
150 h~^ comoving kpc on a side. The dark matter particles in the highest resolution grid 
are 14.48 h~^ Mq and the spatial resolution of the highest resolution grid is 586 h~^ parsecs 
(comoving). Previous work shows that this particle mass resolution is adequate to fully 
resolve the collapse of the halo (Abel et al. 2002; O'Shea & Norman 2007). 

All simulations are performed using the adaptive mesh cosmology code Enzo, which is 
described in Section 12.11 The simulations are started at z = 99 and allowed to evolve until 
the collapse of the gas within the center of the most massive halo, assuming the presence 
of an unevolving soft UV background with intensities in the Lyman- Werner band of Jlw = 
0, 10-24, 10-23■^ 10-22-5, 10-22, 10-21-67^ 10-21-33, and 10-2^ erg s-1 cm-2 Hz-1 sr-1. This 
range covers a much wider range of parameter space than the results described in Machacek 
et al. (2001), and completely encompasses the range of possible Jlw values suggested by 
Wise & Abel (2005), for a wide range of mean Population HI stellar masses. Note that 
many publications use Flw instead of Jlw- Flw has units of erg s-^ cm-2 g^-^, and thus 
Flw = 4:71 Jlw- We will use Jlw throughout this paper, and for convenience express values 
in units of J21, where J21 is the mean intensity of the UV background in the LW band in 
units of 10-21 erg s'^ cm~^ Hz^^ sr^^. 

The equations of hydrodynamics are solved using the PPM method with a dual energy 
formulation, which is required to adequately resolve the thermal properties of gas in high- 
Mach flows. The nonequilibrium chemical evolution and optically thin radiative cooling 
of the primordial gas is modeled as described in Section 12.11 following 9 separate species 
including molecular hydrogen (but excluding deuterium), with an initial electron fraction of 
2.35 X 10-"^ (which is roughly consistent with Peebles (1968) for an = 0.3, Qb = 0.04 
universe). Note that the initial electron fraction in the calculation is relatively unimportant 
to the molecular hydrogen formation rates in halo cores, as the electron fraction at the center 
of a given halo is controlled primarily by mergers and the shock formed by accretion of gas 
onto the halo. 

Adaptive mesh refinement is used throughout the innermost high resolution region such 
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that cells are refined by factors of two along each axis, with a maximum of 22 total levels of 
refinement. This corresponds to a maximum spatial resolution of 115 h~^ AU (comoving) at 
the finest level of resolution, with an overall spatial dynamical range of 5.37 x 10*. To avoid 
effects due to the finite size of the dark matter particles, the dark matter density is smoothed 
on a comoving scale of ~ 0.5 pc (which corresponds to ~ 0.03 proper pc at 2; ~ 18). This is 
reasonable because at that scale in all of our calculations the gravitational potential in the 
halo of interest is completely dominated by the baryons. 

Grid cells are adaptively refined based upon several criteria. We refine on baryon and 
dark matter overdensities in cells of 4.0 and 8.0, respectively. This corresponds to a maximum 
mass of gas or dark matter per cell (on the most highly refined static grid) of M^ax = 12.78 
and 166.16 Mq, respectively. In addition, the MinimumMassForRefinementLevelExponent 
parameter is set to —0.2 for both the dark matter and baryon overdensity refinement criteria, 
meaning that the mass required to refine to a higher level decreases as a function of increasing 
level, as: 

MrefiL) = M^a. X 2-°^^ (1) 

where L is the current level of refinement. The negative exponent means that the 
mass resolution in the calculations is super-Lagrangian - for example, Mref{L = 20) = 
0.0625 Mref{L = 0). In addition to refining on baryon and dark matter overdensity, these 
simulations include refinement criteria which ensure that shocks are always well-resolved, 
that the cooling time in a given cell is always longer than the sound crossing time of that 
cell, and that the Jeans length is always resolved by at least 16 cells. This last criterion 
guarantees that the Truelove criterion (Truelove et al. 1997) is always resolved by a factor 
of four more cells in each dimension than is strictly necessary, ensuring that no artificial 
fragmentation will take place. 

3. Mean halo properties at collapse 

Figures [1] and [2] show projections of baryon density and temperature for the J21 = and 
J21 = 1 calculations, respectively, at the epoch at which each calculation collapses, defined 
as the redshift at which the baryon number density reaches ^ lO^'' cm~'^. Note that due 
to the rapid evolution of gas at high density, the "collapse redshift" depends very weakly 
on the exact choice of density threshold. The calculations are started from the same set of 
initial conditions, and the J21 = 1 calculation is clearly a later stage in the evolution of the 
J 21 = run. The satellite halos surrounding the primary halo which collapses at z=24.12 
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in the J21 = run have merged by z=17.32 when the halo in the J21 = 1 case collapses 
(cf. Fig. [8]). With virial masses of 5.68 X 1O^M0 and 1.26 x lO^Mg, respectively, the latter 
halo is more than 20 times as massive as the former. Though quite different in mass, both 
halos exhibit similar morphologies - they are extremely centrally-concentrated, and only a 
single condensed object (that is to say, a primordial protostellar cloud core) is visible in the 
highest-resolution panel in each image. The existence of a single cloud core is common to 
all calculations discussed in this paper. 
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Fig. 1. — Projections of baryon density and temperature at the final output of the J21 = 
simulation {z = 24.119; = 5.68 x 10^ Mq). Top row: Projected log baryon density. 
Bottom row: Projected, mass-weighted log baryon temperature. Left column: Region 2.13 
kpc (proper) across and deep. Middle column: Region 133.3 pc (proper) across and deep. 
Right column: Region 0.130 pc (proper) across and deep. The middle and right column are 
zoomed images focusing on the high baryon density core which forms in the center of the 
halo. 
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Fig. 2. — Projections of baryon density and temperature at the final output of the J21 = 1 
simulation {z = 17.322; M^r = 1-260 x 10^ M0). Top row: Projected log baryon density. 
Bottom row: Projected, mass-weighted log baryon temperature. Left column: Region 2.13 
kpc (proper) across and deep. Middle column: Region 533 pc (proper) across and deep. 
Right column: Region 0.260 pc (proper) across and deep. The middle and right column are 
zoomed images focusing on the high baryon density core which forms in the center of the 
halo. 
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Figure [3] shows several mean halo quantities for all simulations discussed in this paper, 
including halo collapse redshift as a function of Jlw, virial mass as a function of Jlw: virial 
temperature as a function of Jlw^ and the virial mass as a function of collapse redshift. 
Simulations where the soft UV background is turned on are shown by solid squares, and the 
"control" J21 = calculation is shown by an open square. In plots where Jlw is shown on 
the X-axis, the J21 = simulation is placed at logio Jlvf = —24.5. 

Figure [3] shows that there is a clear relationship between the LW intensity and the 
collapse redshift and virial mass of the halo. A larger LW intensity results in a later collapse 
time and larger virial mass because the halo must be hotter in order to have a cooling time 
which is less than a Hubble time despite the depressed H2 abundance in the high-density gas 
at the halo's center. This is discussed further in Section H] The final mass of the halos in the 
simulations with J21 = 1 is approximately a factor of 20 higher than that in the "control" 
simulation. This is in qualitative agreement with Machacek et al. (2001), who suggest that 
there is a "minimum halo mass" as a function of the strength of the LW background of the 
form 



Mth(M0) = 1.25x10^ + 8.7x10^ (iF^j ^Plw<W~^' (2) 

This threshold mass is plotted in panel (c) of Figure [31 Eq. [2] is only strictly valid over 
the range < Flw = '^t^Jlw < 10^^^, because this was the range simulated by Machacek 
et al. (2001). We see that our points for J21 < 0.1 parallel the threshold curve but at a mass 
approximately four times higher. This difference can be ascribed to two factors. First, the 
threshold mass is the lowest possible halo mass that can collapse, derived from a statistical 
sample. Since we study only one realization which focuses on the most massive halo in the 
box, its mass is bound to be higher than the statistical minimum. Second, Machacek et al. 
(2001) 's criterion for cooling catches halos at an earlier stage of evolution compared to our 
data points, which give the halos' virial masses at the time of central baryon collapse (the 
"collapse redshift"). 

Examination of panels (b) and (c) of Figure [3] show that there is some sort of "phase 
change" between J21 = 10~^'^ and 10~^ causing the halo mass at the time of collapse to 
increase steeply and non-monotonically. The collapse redshift steadily decreases as a function 
of UV increasing background strength (i.e. collapse of gas at the center of the halo is 
delayed). However, panel (c) shows that the virial mass steadily increases with increasing 
UV background strength until J21 = 10~^'^, at which point there is a jump in mass by more 
than a factor of four, above which the halo mass increases only slightly with increasing UV 
background strength. That something interesting should happen at these LW intensities is 
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consistent with the findings of Yoshida et ah (2003) who found that H2 coohng is strongly 
suppressed at J21 = 0.1. They predicted that star formation would be inhibited since the 
equilibrium H2 fraction is below the critical fraction for cooling derived by Tegmark et al. 
(1997). Contrary to these predictions, we find collapse not only at J 21 = 0.1, but also 
at values as large as J21 = 1. We note that highest value considered by Machacek et al. 
(2001) was J21 = 0.0796, and the highest simulated by Yoshida et al. (2003) was J21 = 0.1. 
Therefore our cases J21 > 0.1 have not been examined before, and certainly not at the 
resolution of our simulations. 

How is the gas in the center of these halos able to cool and collapse in such high radiation 
backgrounds? This is analyzed in some detail in the next section. A hint is provided in (d) 
of Figure [31, which plots the halo virial temperature as a function of Jlw- At background 
strengths above J21 = 10~^'^, the virial temperature is consistently approximately 10^ K, 
which is roughly the temperature at which atomic hydrogen cooling is effective. However, an 
examination of the radial temperature profiles shows that high-density gas in the center of 
the halo never reaches these temperatures, but is more typically 2000 K. Atomic line cooling 
is unimportant at such low temperatures; however the H2 cooling rate per particle is roughly 
100 times as large at 2000K as at 500K. H2 cooling still operates in the halo centers despite 
low H2 abundance due to the higher cooling rates and long evolutionary timescales. In the 
rest of this section we merely present additional mean properties of the halo and its central 
region at collapse, and defer discussion of the relevant timescales to Section HI 
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Fig. 3. — Mean halo quantities for several simulations with the same cosmic realization but 
a range of Lyman- Werner molecular hydrogen photodissociating flux backgrounds. Panel 
(a): vs. halo collapse redshift. Panel (b): halo virial mass vs. halo collapse redshift 
Panel (c): halo virial mass vs. Jlw Panel (d): halo virial temperature vs. Jlvf The J21 = 
"control" result are shown as an open square (and is at log J^vf = —24.5 in the panels which 
are a function of ^lw)- In the bottom left panel, the dashed line corresponds to the fitting 
function for threshold mass from Machacek et al. (2001), Eqtn. 8. 
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Figure m shows several properties of the halo core at the epoch of collapse, including the 
core temperature, molecular hydrogen fraction, and spherically-averaged accretion rate as a 
function of ultraviolet background strength, and the spherically-averaged accretion rate as a 
function of the molecular hydrogen fraction. All "core" values are spherically-averaged and 
measured at the mass shell where 100 Mq of gas is enclosed. Panel (a) shows that the core 
H2 fraction decreases significantly with increasing FUV flux, with a corresponding increase 
in the core temperature (Panel (b)). This relationship is similar to that noted in O'Shea 
& Norman (2007), where the amount of molecular hydrogen at densities of 10^ — 10® cm^'^ 
varies between simulations, and correlates strongly with baryon temperature. Given that the 
accretion of gas onto the protostellar cloud is subsonic, this results in a strong relationship 
between the soft UV background flux and accretion rate onto the protostellar cloud, as shown 
by Panel (c), with accretion rates varying by more than a factor of 30 between the J21 = 
and 1 cases. Panel (d) shows the strength of the correlation between core H2 fraction and 
accretion rate. Note that the values discussed above are insensitive to the exact deflnition 
of the halo "core." 
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Fig. 4. — Mean quantities within the central IOOMq core for several simulations with the 
same cosmic realization but a range of LW background intensities. Panel (a): Jlw vs. 
baryon core temperature. Panel (b): Jlw vs. baryon core H2 fraction. Panel (c): Jlw vs. 
instantaneous accretion rate. Panel (d): Baryon core H2 fraction vs. accretion rate. The 
J21 = "control" result are shown as an open square (and is at log Jlw = —24.5 in the 
panels which are a function of Jlw)- All "core" values are spherically-averaged and measured 
at the mass shell where 100 M© of gas is enclosed. 



-17- 



Figure [5] shows the halo core mass and mass fraction at the epoch of collapse (note 
that the definition of "core" is somewhat different than in the previous figure). The halo 
"core mass" is defined as being all gas at a density of nj:^ = 10^ cm~^ or above, and the 
core mass fraction is defined as the core mass divided by the virial mass of the halo at 
the epoch of collapse. This choice of minimum density ensures that we have captured the 
entirety of the quasistatically contracting analog of a galactic molecular cloud core, and for 
practical purposes is comparable to the gas fraction of cold, dense gas discussed in Machacek 
et al. (2001). The size of the halo core increases in absolute terms as the strength of the 
photodissociating background is increased, from ^ 10^ Mq until J21 = 0.1, where it reaches 
~ 11,000 Mq. Above this value of J21, the core mass then begins to decrease, reaching 
^ 7600 M0 at J21 = 1. The core mass fraction generally decreases (though with some noise 
in the relationship) from ~ 1.7 x 10~^ at J21 = to ^ 6.96 x 10""^ at J21 = 1, though the 
total overall change is less than a factor of 2. 
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Fig. 5. — Halo core mass (Panel (a)) and core mass fraction (Panel (b)) for several simula- 
tions with the same cosmic realization but a range of LW background intensities. The "core 



mass" is defined as the mass of gas in the halo which has a density of uh — 10 cm""* or 
above. The "core mass fraction" is defined as the "core mass" divided by the virial mass of 
the halo at the epoch of collapse. The J lw = "control" result arc shown as an open square 
(and is at log Jlw — —24.5 in the panels which are a function of Jlw)- 
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4. Central evolution prior to collapse 

We now examine in more detail tlie evolution of the conditions in the cloud center leading 
to collapse. In Figure[6]we plot the central density, temperature, entropy S = T jnj^ , and H2 
fraction versus time for the calculations with J21 = 10"'^ (solid line), 10~^ (short-dashed line), 
10^^ (long-dashed line) and (dot-dashed line). In these plots, central values of temperature 
and density are defined as being the values of these quantities in the cell with the highest 
baryon density. It is most instructive to compare the two extreme cases, J21 = 10"'^, which 
very closely resembles the J21 = case, and the J21 = 1 case, which according to Yoshida 
et al. (2003) should not collapse at all. In panel (a) we see the central density increase 
modestly from ~ 2 cm^^ to ~ 20 cm~^ over 40 Myr, and then increase rapidly thereafter, 
reaching 10*^ cm~^ in a scant additional 15 Myr. This evolution is driven by H2 cooling, as 
can be seen from the entropy evolution in panel (c). As expected, runaway cooling occurs 
when the H2 fraction reaches ~ 2 x 10~^ (panel d), consistent with the Tegmark et al. (1997) 
analysis. By contrast the J21 = 1 case requires 120 Myr for the central density to increases 
from ~ 2 cm^'^ to ~ 20 cm^^. As shown in panel (c), the first 50 Myr of this contraction is 
adiabatic, followed by an increase in entropy due to some heating event (mergers). The fact 
that central entropy is not decreasing for 90 Myr is indicative of the fact that H2 cooling is 
unimportant over this interval due to extremely low equilibrium H2 fractions {iH2 < 10~^; 
(panel d)). 

However, something interesting happens in the J21 = 1 calculation at t ~ 170 Myr. 
Prior to that time the central temperature has crept up to 5000 K due to the increased virial 
mass. The cooling rate per H2 molecule is about 3 orders of magnitude higher at 5000K than 
at 500K - a typical temperature at the halo center in low UV background evolutions. We 
believe this high temperature is due to a combination of adiabatic heating as the potential 
well deepens, and merger-induced shock heating. As shown in Figure [H three large halos 
merge between z = 20 and z = 18, quadrupling the halo mass (Figure [3l panel (b)). The 
elevated cooling rates cool the gas at the halo center toT ^ 2000if, allowing the gas density 
to increase slightly. The two effects reduce the central entropy from S ~ 700 K cm^ to ~ 
250 K cm^. This is followed by another heating event at t=190 Myr, presumably due to 
another merger, followed by a second more catastrophic cooling event. At t=200 Myr, the 
halo core collapses, driving the central density higher, temperature and entropy lower, and 
H2 fraction higher. 

Why does this collapse occur? This is analyzed in Figure [71 where we plot the evolution 
of the important timescales in the center of the halo. In the J21 = 10~^ case, the 2-body H2 
formation timescale is always roughly an order of magnitude smaller than the photodissoci- 
ation timescale. Consequently, the H2 abundance is out of equilibrium and grows steadily 
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with time (Figure [6] panel d). The H2 coohng time in the central zone (< 1 Myr) is far less 
than the Hubble time (~ 100 Myr), and before too long runaway collapse occurs. At higher 
values of Jlw: the 2-body H2 formation time always hovers around the photodissociation 
timescale, regardless of how long that timescale is. This tells us that H2 is in equilibrium: the 
2-body formations balance the photo-destructions by the LW background. The equilibrium 
H2 fraction is given by 



^ -J— (3) 

where k^- is the rate coefficient for the formation of H~ (the limiting reaction for 
H2 formation at the temperature and density range considered here), Ue is the electron 
density, and k^iss = 1.1 x lO^F^vy*"^ is the H2 photodissociation rate due to the Solomon 
process (Abel et al. 1997; Yoshida et al. 2003). Using T = 2000 K and = 10-^ we get 
fH2 = 1-46 X 10~^ (1^^) ' which is in good agreement with the H2 fractions we see in 
Figure [61 Since 1^2 = nH2/ ^H2 oc fH2^ the smaller fH2 is, the shorter its formation time. 

Now let us consider the H2 cooling time evolution (short-long-dashed line). It is far lass 
than a Hubble time for J21 < 10^^, ensuring that these halo cores will eventually cool and 
collapse. In the J21 = 1 case, the cooling time drops below the Hubble time for the first 
time at t=170 Myr, and again at t=190 Myr where it begins a steady decrease toward 10^ 
yr. Despite the low equilibrium H2 fraction, if we wait long enough the inexorable press of 
time eventually estabhshes the condition tcooi < tuubUe and we get runaway collapse, even 
for J21 = 1. Note that we get catastrophic cooling because we are in the density regime 
(n < 10^ cm~^) where H2 cooling is proportional to the square of the gas density, not linear 
in the gas density (Abel, Bryan, & Norman 2000; Abel et al. 2002). 
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Fig. 6. — Evolution of the central density (Panel (a)), temperature (Panel (b)), entropy 
(Panel (c)), and H2 fraction (Panel (d))as a function of time for the cases J21 = 10~^ (solid), 
10~^ (short-dashed), 10~^ (long-dashed) and 1 (dot-dashed). 
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Fig. 7. — Evolution of timcscalcs in the central zone for the cases J21 = 10^^, 10^^, 10~^ 
and 1 (Panels (a)-(d), respectively): 2-body H2 formation time (solid); 3-body H2 formation 
time (short-dashed); H2 coUisional dissociate time (long-dashed); H2 photodissociation time 
(dot-dashed); H2 cooling time (short-long-dashed). 
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Fig. 8. — Merger history of the M^r = 1.26 x lO'^M©, T^r = 10^ K halo that collapses at 
z=17.32 with a Lyman- Werner background flux J21 = 1. In comic strip order, the redshifts 
are z=22, 21, 20, 19, 18, and 17.32. Field of view is 53.571 comoving kpc (2.678 proper kpc 
at z = 19). Logarithm of the projected baryon density (column density) is displayed, with 
color table scaled to the maximum and minimum values in each image. 



5. Halo properties at the epoch of collapse 



Figure [9] shows radial profiles of baryon number density, temperature, enclosed mass, 
circular velocity, and RMS Mach number as a function of radius, as well as specific angular 
momentum as a function of enclosed mass, for all simulations discussed in this paper. To 
facilitate comparison, results from the output of each simulation where the peak baryon 
number density at the center of the halo is approximately 10^° cm~^ are shown, in order 
to capture each halo at a similar evolutionary stage, rather than at a fixed point in time. 
Panel (a) shows that all of the simulations have similar density profiles. The scatter in 
number density at a given radius is readily explained by variation in halo mass - halos which 
collapse at later times (due to higher UV background strength) are more massive, and thus 
have higher overall baryon densities. This is shown in another way in Panel (c), in a plot of 
enclosed mass as a function of radius which shows that all halos have very similar profiles 
at r < 10~^ pc, but a significant variation at larger radii which is related to halo mass. 
The plot of temperature as a function of radius in Panel (b) shows an interesting trend - 
as the LW background is increased, the overall halo temperature as well as the halo core 
(where "core" is roughly defined as gas within ~ 1 pc of the halo center) temperature go 
up. In the outskirts of the halo, where baryon densities are low and thus cooling times 
are long, this is due primarily to the increase in halo mass, with the peak temperature 
corresponding approximately to the halo virial temperature. In the halo core, however, this 
temperature is correlated more strongly with the H2 fraction, and thus the UV background, 
as suggested by Figure HI The plots of specific angular momentum as a function of enclosed 
mass and Keplerian velocity fraction (defined as the local circular velocity divided by the 
Keplerian orbital velocity due to mass within that radius) as a function of radius in panels 
(d) and (e) show some mild trends trends. It appears that halos in simulations with stronger 
UV backgrounds tend to have less specific angular momentum at a given mass shell, and 
also tend to have a lower Keplerian velocity fraction at a given radius than gas in halos 
which form in the presence of a lower UV background - indeed, the halo which forms in the 
J21 = 1 calculation has the least angular momentum of all of the simulations, and the lowest 
Keplerian velocity fraction out to ~ 1 pc. We speculate that this is due to transport and 
segregation of angular momentum by turbulence - halos whose collapse is delayed have more 
time for these processes to act within the halo core, resulting in less angular momentum 
overall when the halo ultimately collapses (note that this will be discussed more fully in 
a later paper). The plot of RMS Mach number as a function of radius in panel (f) also 
shows no obvious trend, though at radii below ~ 10~^ pc, simulations with a smaller UV 
background generally have higher RMS Mach numbers. This is due primarily to the gas 
being colder, and not the gas velocities being higher. 
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Fig. 9. — Evolution of several spherically-averaged baryon quantities for simulations with 
all values of J2i- Panel (a): p^r^ as a function of radius. Panel (b): Temperature as a 
function of radius. Panel (c): Enclosed gas mass as a function of radius. Panel (d): Specific 
angular momentum as a function of enclosed mass. Panel (e): Keplerian velocity fraction as a 
function of radius. Panel (f): Gas RMS Mach number as a function of radius. All quantities 
except enclosed gas mass are mass-weighted, and all simulations are shown at the point 
where the maximum number density is approximately 10^° cm~^. Line types and weights 
correspond to simulations, as follows. Black solid line: J21 = 0. Black short-dashed line: 
J21 = 10"^. Black long-dashed line: J21 = 10~^-^. Black dot short-dashed line: J21 = 10"^. 
Black dot long-dashed line: J21 = 10^^'^. Black short dashed-long dashed line: J21 = 10^^. 
Red solid line: J21 = 10"°-^''. Red short-dashed line: J21 = 10~°'^'^. Red long-dashed line: 
J21 = 1. 



- 26 - 



Figure [TO] shows the H2 fraction, electron fraction, H~ fraction, ratio of coohng time to 
sound crossing time, ratio of coohng time to dynamical time, and ratio of sound crossing 
time to dynamical time as a function of radius, for all simulations. Outputs and line types 
correspond to those in Figure O There is a clear relationship at all radii between H2 fraction 
and UV background strength in panel (a) - as the FUV background is increased, the overall 
H2 fraction decreases. This difference is most noticeable at radii of 0.1 — 10 pc, but is 
maintained at larger and smaller radii. There is a "kink" at approximately 10^^ pc where 
the values are all quite similar. This corresponds to a baryon number density of ~ 10*^ cm~'^ 
in all simulations, which is where 3-body H2 formation begins to occur. Panels (b) and (c) 
show that e~ and H~ fractions track each other, which is to be expected - the local electron 
fraction controls the amount of H~ which can be produced, which is the limiting reaction in 
the dominant mode of H2 formation for n < 10^ cm~^ - and there is a general trend with 
increasing electron and H~ fractions with increasing UV background flux. 

The plots of the ratio of cooling time to sound crossing time, ratio of cooling time to 
dynamical time, and ratio of sound crossing time to dynamical time as a function of radius 
shown in panels (d), (e) and (f) display some interesting trends. The three simulations 
with the lowest UV background strengths have a much longer cooling time than the other 
calculations, which are all grouped roughly together with no discernible pattern. This agrees 
well with the plot of the ratio of cooling time to dynamical time, where a similar trend is 
observed. The difference between the three simulations with the lowest UV background 
strengths and the others is due to the somewhat higher H2 fraction in the halo core of 
these calculations. The gas temperature is ~ 200 K at the center of the halo cores in these 
calculations, and the gas cannot cool further, which results in a very long cooling time. 
The sound crossing time is also increased, but not as significantly. All simulations have 
approximately the same ratio of sound crossing to dynamical times at all radii. Given that 
the cooling time is typically longer than both the sound crossing time and dynamical time, 
one can infer that the collapse of the halo is occurring quasi-statically for all simulations. 
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Fig. 10. — Several spherically-averaged baryon quantities for simulations with all values of 
J21. Panel (a): Molecular hydrogen fraction as a function of radius. Panel (b): Electron 
fraction as a function of radius. Panel (c): H~ fraction as a function of radius. Panel (d): 
Ratio of gas cooling time to sound crossing time as a function of radius. Panel (e): Ratio 
of gas cooling time to dynamical time as a function of radius. Panel (f): Ratio of gas sound 
crossing time to dynamical time as a function of radius. All quantities except enclosed gas 
mass are mass-weighted, and all simulations are shown at the point where the maximum 
number density is approximately 10^° cm~^. Line types and weights correspond to those in 
Figure O 
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Figure [TT] shows radial velocity as a function of radius, accretion time (defined as 
Mfy^enc/M) as a function of enclosed baryon mass M as a function of enclosed baryon mass 
and M as a function of time for all simulations. Output time and lines correspond to those 
in Figures [3 and Uni Panel (a) shows a clear relationship between the strength of the soft UV 
background and the radial velocity, with the calculations that have higher UV backgrounds 
also having higher infall velocities. This is easily understood by the quasistatic contraction 
of the gas, which takes place at or below the local sound speed in the gas, as discussed by 
Abel et al. (2002) and O'Shea & Norman (2007). Given that the halo temperatures are 
systematically higher in calculations where the UV background is stronger, the sound speed 
is higher, and thus the rate at which the gas at the center of the halo contracts is higher. 
This leads to lower accretion times (Panel (b)) and higher accretion rates (Panel (c)), with 
the accretion rates for the simulations with the strongest UV backgrounds being higher than 
those with a low (or no) UV background. Panel (d) shows that the time evolution of the 
accretion rate in all simulations is qualitatively similar, but varies strongly in absolute mag- 
nitude and in the time at which the accretion rate onto the evolving protostellar cloud peaks, 
suggesting that the stars forming in these simulations may have very different evolutionary 
histories. 
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Fig. 11. — Several spherically-averaged baryon quantities for simulations with all values of 
J21. Panel (a): Radial velocity as a function of radius. Panel (b): Accretion time as a function 
of enclosed baryon mass. Panel (c): Instantaneous accretion rate as a function of enclosed 
baryon mass. Panel (d): Estimated accretion rate as a function of time. All quantities except 
enclosed gas mass are mass-weighted, and all simulations are shown at the point where the 
maximum number density is approximately 10^° cm~^. Line types and weights correspond 
to those in Figure [91 The upper and lower light short-long-dashed curves which extend from 
the upper left corner of Panel (b) correspond to the main sequence lifetime of a massive 
Population III star of that mass and the Kelvin-Helmholtz timescale of a Population III 
timescale with a given luminosity and radius. All values are taken from Schaerer (2002). The 
three light diagonal short-dashed lines which extend from bottom left to top right in panel 
(b) correspond to masses accreted using constant accretion rates of (from top to bottom) 
m = 10~^, 10~^, and 10~^ M0/yr. The light horizontal short-dashed line in panels (c) and (d) 
correspond to the "critical" accretion rate of Omukai & Palla (2003), m ~ 4 x 10~^ Mq/jt. 
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6. Halo properties at fixed redshift 

Figure [12] sliows tlie status of several spherically-averaged quantities for halos from all 
simulations at a fixed redshift, z = 25, which is shortly before the baryons in the most massive 
halo in the J21 = simulation collapses to high density. This is at an early point in the 
evolution of the majority of these halos, and several conclusions can be drawn from the plots 
of number density, temperature, H2 fraction, and radial velocity as a function of radius. The 
molecular hydrogen fraction declines monotonically with increasing UV background strength, 
in agreement with an estimate for equilibrium H2 fraction at low {n <^ 10® cm~^) densities: 



JH2 - -r (4) 

where k^- is the rate coefficient for the formation of H~ (the limiting reaction for 
H2 formation at the temperature and density range considered here) and k^w is the rate 
coefficient for H2 in the LW band (Abel et al. 1997; Yoshida et al. 2003). Halo cores at higher 
central densities are departing from the approximations used in the above estimate and thus 
have values of the H2 fraction that are not quite the equilibrium values. Halo cores which 
have reached core H2 fractions higher than ~ 10~^ have been able to cool to significantly 
below the virial temperature (since the cooling time is much less than the Hubble time at 
that value of the H2 fraction), allowing the baryon density to increase. The effect of cooling 
can also be seen in the plot of radial velocity as a function of radius, where the core regions 
of halos that have cooled are contracting (albeit very slowly), while those which are still at 
the virial temperature are not showing significant signs of contraction within ~ 20 pc. All 
halos have some infiow at large radii, due to the infall of gas from adjoining filaments. 
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Fig. 12. — Several spherically-averaged baryon quantities for simulations with all values 
of J21. Panel (a): Number density as a function of radius. Panel (b): Temperature as a 
function of radius. Panel (c): H2 fraction as a function of radius. Panel (d): Radial velocity 
as a function of radius. All quantities except enclosed gas mass are mass-weighted, and all 
simulations are shown at z = 25 (shortly before the J21 = simulation collapses to high 
density. Line types and weights correspond to simulations, as follows. Black solid line: 
J21 = 0.0. Black short-dashed line: J21 = lO^'^. Black long-dashed line: J21 = 10^^'^. Black 
dot short-dashed line: J21 = 10^^. Black dot long-dashed line: J21 = 10~^'^. Black short 
dashed- long dashed line: J21 = 10~^. Red solid line: J21 = 10"°'^^. Red short-dashed line: 
J21 = 10~°-^^. Red long-dashed line: J21 = 1. 
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7. Evolution of two representative halos 

In this section we compare tlie evolution of two representative simulations. We choose 
the calculations which correspond to J 21 = and 1. These two simulations represent the 
extremes in halo evolution - the J21 = calculation is an example of the evolution of a 
"standard" halo, as discussed in previous literature, and the J21 = 1 simulation is the most 
extreme example of a Population III protostar forming in a Tyir ~ 10^ K halo available from 
our suite of calculations. As described in Section 12.21 both calculations start from the same 
initial conditions, with only the strength of the LW background being different. 

Projections of the log baryon density and temperature for the J21 = and 1 simulations 
are shown in Figures [1] and [2l respectively. These projections are at the epoch of collapse, 
which is z = 24.1 (17.3) for the J 21 = (1) simulation. The halo masses are significantly 
different: 5.68 x 10^ and 1.26 x 10'' Mq for the two calculations, refiecting the rapid pace of 
mergers at that redshift, and the amount of structure evolution that takes place over that 
relatively short time period (~ 8.2 x 10^ years). This can clearly be seen in the left column 
of both figures, when the large amount of structure apparent in Figure [1] has merged into 
the main halo by the time of protostellar cloud formation in the J21 = 1 case, as seen in 
Figure [21 

Despite the differences in evolutionary states, the halos themselves are quite similar in 
Figures [1] and [2l Both calculations show halos that are strongly centrally concentrated, as 
shown by the center column in both figures. This column shows projections of the halos 
which are centered on the baryon peak and are approximately one virial radius across (and 
thus are scaled differently in both images). No fragmentation of the halo core is visible. The 
rightmost panel in both calculations shows the center of the halo, where the protostellar 
cloud is evolving. There is still no evidence for fragmentation up to a central baryon density 
of n ~ 10^° cm~^. 

Figures [T3] and [H] show the temporal evolution of several spherically-averaged baryon 
quantities for the J21 = and 1 calculations, respectively. These plots show baryon number 
density, enclosed baryon mass, temperature, H2 fraction, ratio of cooling time to sound 
crossing time, and ratio of cooling time to dynamical time, as a function of radius, evolving 
from a number density of ~ 10 cm~^ to ~ 10^° cm~^ in both calculations. Some qualitative 
commonalities are obvious between the two calculations. Both halos steadily grow in central 
density (Panel (a)), which is coupled to the growing molecular hydrogen fraction in the 
central region of the halo (Panel (d)). Both halos collapse quasi-statically, as shown in 
panels (e) and (f). Also, the temperature profiles (as shown in Panel (c)) are broadly 
similar, in that the gas in the halo central region cools to some minimum value at r ~ 1 — 2 
pc (corresponding to n ~ 10'^ cm~^ in both simulations) and then creeps steadily upward as 
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gas collapses to higher densities. 

Though there are quahtative similarities between the evolution of the two calculations, 
detailed examination shows some significant quantitative differences. The time that the halo 
in the J21 = simulation takes to evolve from a central density of n ~ 10 cm~^ to 10^° cm~^ 
is ~ 24 Myr, with 21 Myr of that being the time required to evolve from n ~ 10 to n 

~ 10'* cm"^. The time required for the J21 = 1 calculation to evolve from n ~ 10 cm~^ to 
lO"*^*^ cm^"^ is slightly shorter, ^ 17 Myrs, with 8 Myrs required to evolve to n ^ 10'* cm^"^. The 
temperature evolution is also somewhat different - the J21 = simulation has a minimum 
temperature at r ~ 1 pc of 200 K, with n ~ 10*^'^ and a molecular hydrogen fraction of 
~ 10^^, and roughly 1000 Mq of enclosed gas. As the gas within this radius evolves to 
higher densities, the temperature increases by a factor of a few, peaking at T ~ 500 K at 
n ~ 10^*^ cm~^, with in^ — 0.02. The gas coohng time within r ~ 1 pc in the J21 = 
simulation is always larger than the sound crossing and dynamical times by at least a factor 
of two, implying a quasistatic gas contraction at all times. 

The J21 = 1 simulation evolves somewhat differently. The collapsing gas reaches a 
minimum temperature of ~ 800 K at r ~ 2 pc, with a baryon density of ~ 2 x 10^ cm~^, 
an enclosed gas mass of roughly 10* M©, and a molecular hydrogen fraction of ~ 10~^, two 
orders of magnitude below that in the J21 = calculation at the equivalent temperature 
minimum. The low H2 fraction is due to the high LW background radiation flux. As the 
halo evolves to higher densities, the gas temperature and H2 fraction also creep upwards. 
The temperature of the gas when the peak reaches roughly 10*° cm~^ is ~ 1200 K, with a 
H2 fraction of ~ 4 x 10"'^. The cooling time at r ~ 2 pc is an order of magnitude higher than 
the sound crossing or dynamical times, implying an extremely slow contraction of the gas at 
that radius. As the central density increases, however, the cooling time becomes lower, with 
T^cooi/T^dyn ~ 0.7 and Tcooi/T^cross dipping to ~ 0.8 at its minimum, but the latter increasing 
toward the central density peak to a ratio of ~ 10. This implies a more rapid contraction 
of gas in the halo center in the simulation with a higher UV background, despite the lack of 
molecular hydrogen. It is also worth noticing that the J21 = 1 calculation has H2 fractions 
at n ~ 10'* and 10*^ cm~^ of ~ 2 x 10^^ and 8 x 10^*, respectively, as compared to f/j-j — 10~^ 
and 2 x 10^"^ at the same densities in the J21 = calculation. 
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Fig. 13. — Evolution of several spherically-averaged baryon quantities as a function of time 
for the J21 = simulation. Panel (a): Number density as a function of radius. Panel 
(b): Enclosed gas mass as a function of radius. Panel (c): Temperature as a function of 
radius. Panel (d): H2 fraction as a function of radius. Panel (e): Ratio of cooling time to 
sound crossing time as a function of radius. Panel (f): Ratio of cooling time to dynamical 
time as a function of radius. All quantities except enclosed gas mass are mass-weighted. 
Line types and weights correspond to different times, as follows. Red short dashed-long 
dashed line: z = 27.678, t = 1.1088 x 10^ years. Red dot long-dashed line: z = 26.010, 
At = 1.0428 X 10^ years. Red dot short-dashed line: z = 24.565, At = 1.0428 x 10^ 
years. Red long-dashed line: z = 24.273, At = 2.291 x 10^ years. Red short-dashed line: 
z = 24.171, At = 8.0959 x 10^ years. Red solid line: z = 24.148101, At = 1.9034 x 10^ 
years. Black short dashed-long dashed line: z = 24.132556, At = 1.2528 x 10^ years. Black 
dot long-dashed line: z = 24.126839, At = 4.6131 x 10^ years. Black dot short-dashed 
line: z = 24.124178, At = 2.1483 x 10^ years. Black long-dashed line: z = 24.121760, 
At = 1.9517 X 10^ years. Black short-dashed line: z = 24.120224, At = 1.2404 x 10^ years. 
Black solid line: z = 24.119543, At = 5.5024 x 10^ years. 
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Fig. 14. — Evolution of several spherically-averaged baryon quantities as a function of time 
for the J21 = 1 simulation. Panel (a): Number density as a function of radius. Panel (b): 
Enclosed gas mass as a function of radius. Panel (c): Temperature as a function of radius. 
Panel (d): H2 fraction as a function of radius. Panel (e): Ratio of cooling time to sound 
crossing time as a function of radius. Panel (f): Ratio of cooling time to dynamical time 
as a function of radius. All quantities except enclosed gas mass are mass-weighted. Line 
types and weights correspond to different times, as follows. Red short dashed-long dashed 
line: z = 18.384605, t = 1.995 x 10^ years. Red dot long-dashed line: z = 17.737221, 
At = 5.2141 X 10^ years. Red dot short-dashed line: z = 17.394373, At = 1.7253 x 10*^ 
years. Red long-dashed line: z = 17.373929, At = 3.6030 x 10^ years. Red short-dashed line: 
z = 17.330286, At = 7.7255 x 10^ years. Red solid line: z = 17.325077, At = 9.25088 x 10^ 
years. Black short dashed-long dashed line: z = 17.323130, At = 3.4610 x 10^ years. Black 
dot long-dashed line: z = 17.322398, At = 1.3001 x 10^ years. Black dot short-dashed 
line: z = 17.322050, At = 6.1828 x 10^ years. Black long-dashed line: z = 17.321880, 
At = 3.0262 X 10^ years. Black short-dashed line: z = 17.321796, At = 1.4813 x 10^ years. 
Black solid line: z = 17.321757, At = 7.0496 x 10^ years. 



-36- 



8. Neglected Physics and Possible Numerical issues 

In this paper we have examined aspects of the formation of Population III stars in the 
presence of a soft ultraviolet background, using an adaptive mesh refinement cosmological 
structure formation code. Given the nature of the tool and current limitations in computing 
power, some physics were neglected, and possible numerical issues may arise. We discuss 
this here. 

The primordial chemistry model used in these calculations ignores the effects of deu- 
terium, lithium, and the various molecules that form between these elements and ordinary 
hydrogen. Deuterium and lithium have been shown to be unimportant in the temperature 
and density regimes that we have examined in this paper (Galli & Palla 1998, 2002; Lipovka, 
Nunez- Lopez, & Avila- Reese 2005). However, it is possible that they may be relevant in other 
situations of importance to Population III star formation - in particular, regions which have 
been ionized to very high electron fractions may experience significant cooling from the HD 
molecule, which due to its permanent dipole moment makes it more than 100 times more 
effective as a cooling agent than molecular hydrogen (per molecule), and has the potential 
to cool gas down to approximately the temperature of the cosmic microwave background, 
which scales with redshift as Tcmbiz) = 2.73 x {1 + z) K (Flower et al. 2000; GaUi & Palla 
2002; Lipovka et al. 2005). This gives a minimum baryon temperature of approximately 55 
Kelvin at 2; = 20 and could reduce the minimum accretion rate onto a primordial protostar 
significantly. Lithium, while in principle an effective coolant as LiH, is safely ignored, since 
only a tiny fraction of lithium is converted into LiH (Mizusawa, Omukai, & Nishi 2005). A 
final chemical process that is omitted from our calculation is heating caused by molecular 
hydrogen formation at high densities (n > 10^ cm"'^), the inclusion of which may result in 
differences in the temperature evolution of the gas in the highest density regimes consid- 
ered here. This does not significantly affect the conclusions of our paper, as the significant 
differences we observe between simulations occur at densities much lower than ~ 10^ cm~^. 

Self-shielding of the photodissociating background by molecular hydrogen in the high- 
density gas is ignored in these calculations. Though this effect could in principle be impor- 
tant, the actual column densities of molecular hydrogen are typically far too small to actually 
block the soft UV flux. According to Machacek et al. (2001) (and references therein), a col- 
umn density of 5 x 10^^ cm~^ is enough for shielding to become important. However, this was 
derived for a static distribution of H2, while the LW band consists of hundreds of individual 
lines whose width in this case is dominated by Doppler broadening. It is useful to note in this 
case that the average line width is ~ 2 km/s and the RMS baryon velocity in our calculations 
are ~ 4 km/s. In order for self-shielding to be important in the case of a turbulent medium, 
the column density must be much higher. Typical maximum H2 column densities in our cal- 
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culations are on the order of 10^^ cm~^, but these occur late in the collapse of the halo core, 
and in the highest density regions the cooling and H2 production times are much shorter 
than the photodissociation time scale, at which point self-shielding becomes unimportant. 
It is worth noting that there are some regimes where self-shileding can be critical. Susa 
(2007) finds that self-shielding can strongly affect the evolution of a collapsing primordial 
halo which is being illuminated with photodissociating flux by a neighboring star, though 
the situation is somewhat idealized. Additionally, we do not consider the more complicated 
effects relating to Population III stars which form in halos in the cosmic neighborhoods 
where previous generations of stars have existed. This allows us to ignore complex radiative, 
chemical and dynamical effects that would vastly complicate our calculations (O'Shea et al. 
2005a; Johnson, Greif, & Bromm 2006; Yoshida et al. 2006a; Greif et al. 2007; Abel, Wise, 
& Bryan 2007; Ahn & Shapiro 2007). 

A further effect that is ignored in this paper is H~ photodetachment. This could in 
principle be a significant effect - H~ is the catalyst for molecular hydrogen formation in the 
dominant H2 formation channel at densities < 10® cm~^, and it can be detached by photons 
with energies hv > 0.75 eV. Photons with this energy would be produced in great numbers by 
the same massive Population III stars that we assume are producing the molecular hydrogen 
photodissociating background. However, as shown by Machacek et al. (2001), the rate of 
photodetachment is orders of magnitude smaller than the rate of H~ formation at densities 
comparable to that found in the centers of the halos examined in this work, and thus the 
process of H~ photodetachment can be safely ignored. 

In this paper, we examine the effects of molecular hydrogen dissociating backgrounds 
which range in strength from J^vf = 10~^^ to 10~^^ ergs~^ cm~^ Hz~^ sr~^ in the LW band 
(11.18 — 13.6 eV). This range is in good agreement with the expected range of photodisso- 
ciating backgrounds predicted by Wise & Abel (2005), but it is still finite. Examination of 
our simulations show that the case with the lowest UV background (J21 = 10^^) is almost 
identical to the "control" (J21 = 0) calculation, justifying our choice of minimum value. 
Our upper value is consistent with Wise & Abel, and examination of Figures [3] and [9HTT] 
suggests that increases in the strength of the UV background (within reasonable values) will 
result in a further delay in halo collapse, but no major change in the mode of star formation 
observed. This is an effect of the halo properties - a small amount of molecular hydrogen 
will always exist in these halos, allowing the gas at the center to cool and contract quasi- 
statically to higher densities. Once the gas has collapsed to very high densities, rapid II2 
formation via the 3-body process will occur, essentially independent of the strength of the 
soft UV background, and the gas will be able to cool down to ~ 200 K very quickly. 

The effects of magnetic fields are completely ignored in the simulations discussed in this 
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work. Magnetic fields are discussed in detail in Paper I, but, to summarize, a fairly high 
seed magnetic field is needed to be dynamically significant at the relatively low densities we 
explore in this work. The possible importance of magnetic fields has been explored in analytic 
and semi-analytic work (Tan & McKee 2004; Silk & Langer 2006; Maki & Susa 2007). We 
will examine the possible evolution of magnetic fields within the context of cosmological 
AMR simulations in a later paper (O'Shea & Turk 2007, in preparation). 

The simulations presented in this paper are generated assuming a cosmology that is 
somewhat different than the currently-favored WMAP Year 111 "best-fit" model (Spergel 
et al. 2006). Most importantly, our value of ag is 0.9, which is significantly higher than the 
WMAP value of 0.761. The general effect of a higher as is to cause structure formation to 
take place earlier, and thus at a given redshift one would expect significantly more halos in 
our chosen cosmological model than in the WMAP Year III model. However, the evolution 
of any single cosmological halo, such as the one examined in this work, is not particularly 
affected by this parameter in the sense that we are not examining halo stastical properties. 
In addition to ag, the ratio of Qb/^m in our simulations is 0.1337. while it is 0.1746 in the 
WMAP Year III cosmology. This may result in some small quantitative differences in, e.g., 
the redshift of halo collapse, but should not significantly affect our results. 

We direct the reader to Paper I for a detailed discussion of other possible numerical 
issues, such as the generation of cosmological initial conditions, the assumption that the 
cold dark matter model is correct, the choice of halo in our simulations, and the size of the 
simulation volumes used. 



9. Discussion 

This paper explores the formation of Population III stars in simulations with a constant 
soft UV background. Our results agree well qualitatively with that of Machacck ct al. (2001); 
we both find that a soft UV background can delay the formation of molecular hydogen, and 
thus the cooling and collapse, of small (~ 10^ M©) cosmological halos in which Population III 
stars form. We also find that increasing the soft UV background increases the minimum halo 
mass required for a halo to collapse (in a way similar to that of smoothing the dark matter 
power spectrum at small scales - see O'Shea & Norman (2006) for a discussion). Machacek 
et al. derived a mass threshold for collapse as a function of the LW background flux that 
agrees well with our simulations, though the halo masses in our calculation are significantly 
higher. This is due to the halo that we examine being a "typical" halo rather than at the 
threshold mass for star formation. Presumably, if we performed these calculations using 
many halos in a range of cosmological realizations, we would find a minimum halo mass 
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that is somewhat lower than the masses seen in our calculations. Our work is a significant 
improvement upon that of Machacek et al. in some respects, as our simulations are much 
more highly resolved and we examine the evolution of a single halo over a much wider range 
of soft UV background fluxes. One drawback of our work compared to Machacek et al. is 
that we examine the evolution of a single halo, albeit with a broad range of UV backgrounds. 

Our work is similar to the calculations in Yoshida et al. (2003) which examine Population 
111 star formation in the presence of a photodissociating ultraviolet background. They flnd 
that for values of Jlw above 10^^"^ erg s^^ cm^^ Hz^^ sr^^, hydrogen molecules are rapidly 
dissociated and gas cooling is inefficient, implying that the collapse of gas in the halo center 
may be delayed until atomic line cooling dominates, in halos with T^j^ ^ 10^ K. While we 
see delays in halo collapse, we do see that a tiny amount of H2 can form even at high FUV 
background strengths, eventually allowing gas to cool and collapse to high densities without 
ever reaching temperatures at which atomic line coofing would be effective (T > 10^ K). 
Yoshida et al. make this statement based on a set of simulations using two different \"alucs 
of Jlw (10~^^ and 10^^^), which bracket the "break" in properties seen in our results. The 
simulations used to obtain this result had SPH particles with m^^s = 100.0 h~^ M© and 
thus an effective mass resolution comparable to the mass of gas within the quasistatically- 
coUapsing central region of a halo. Our mass resolution is more than two orders of magnitude 
higher in the center of the cosmological halos that we examine. The differences between our 
result and theirs is likely due to the significant differences in resolution and may in fact be 
consistent once resolution is taken into account. 

We see a strong relationship between accretion rates onto the protostellar cloud and 
the strength of the photodissociating background. Calculations with higher ultraviolet back- 
ground strengths typically have a larger spherically-averaged baryon accretion rate onto 
the evolving protostellar cloud. This is clearly due to variation in halo central temper- 
ature relating to the amount of molecular hydrogen existing in the center of the halo at 
the epoch of collapse. Accretion onto the evolving protostellar cloud is subsonic, and thus 
regulated by the local sound speed (which scales as T°'^). This implies some relationship 
between the strength of the photodissociating background and the final stellar mass of the 
primordial star - however, the details of this relationship are unclear, and depends on many 
factors. For example, Omukai & Palla (2003) suggest that an increase in accretion rate above 
m ~ 4 X 10~^ Mo/yr will actually result in a decrease in the final mass of the star due to 
radiative feedback from the evolving protostar. However, their results use one-dimenstional, 
fairly idealized models, and geometrical effects may be important. This is explored in Tan & 
McKee (2004), who use a combination of analytic and semi-analytic models of the evolving 
system. They suggest that accretion onto the protostar is highly non-spherical, and is in fact 
mediated by an accretion disk. They argue that a larger accretion rate onto the disk will 
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lead to a larger star, though with some hmitations. Finally, Omukai (2001) and Omukai & 
Yoshii (2003) study the IMF of stars that form in Tyj^. > 10^ K metal-free protogalaxies in 
the presence of an H2 photodissociating background. They conclude that photodissociation 
actually decreases the Jeans mass of the gas at high densities, and thus reduces the fragmen- 
tation mass scale of the clouds and presumably the stellar mass. They use one-dimensional, 
spherically symmetric simulations, however. We can explore this issue in more detail in later 
calculations, when appropriate models for shielding, cooling via H~ , and more advanced 
chemistry have been implemented into Enzo. This will be examined in a later paper. It 
is worth noting that the range of accretion rates observed in this situation are within the 
range of rates seen in the ACDM simulations discussed in Paper I. This imphes that the 
mass ranges of the resulting stars will not be significantly different than in the "standard" 
Population III star formation scenario. 

We see no fragmentation of the high-density baryon core, up to a baryon number den- 
sity of n ~ 10^'' cm~^, for the entire range of simulations explored in this paper. This is 
due to a combination of effects, but predominantly the poor cooling properties of molecular 
hydrogen. The gas is relatively hot (~ 1000 Kelvin) and thus has a high sound speed, which 
helps to damp out perturbations in the halo center which would otherwise result in multi- 
ple fragments. This effect is exacerbated in simulations with high ultraviolet background 
strengths, as the temperatures (and thus the sound speeds) are generally higher. This result 
implies that there is no fundamental change in the mode of Population III star formation as 
halos grow, and that more massive "proto-galactic" halos, with T^j^ ^ 10^ K, will continue 
to form a single massive star per halo. These halos will have a much larger binding energy 
than the smaller, M ~ 10^ — 10^ M© halos which have traditionally been examined by nu- 
merical simulations, and implies that multiple generations of star formation may be able to 
take place in a single halo. This leads to the possibility of "self-enrichment," where a single 
Population 111 star enriches the high-density center of the halo to metallicities high enough 
to change the cooling properties of the gas. This could cause a strong change in the IMF 
(Santoro & ShuU 2006; Smith & Sigurdsson 2007), and these objects could make a significant 
contribution to the reionization history of the universe (Mackey, Bromm, & Hernquist 2003; 
Greif & Bromm 2006). We will examine whether the high-density gas in the centers of these 
more massive halos fragments in a later work. 

The values of accretion rates onto the evolving protostellar cloud which are observed 
with increasing photodissociating background flux (which are always comparable to those 
observed in Paper I), in addition to the apparent lack of fragmentation of gas in the halo 
center up to densities of n/j ~ 10^° cm~^, implies that there is no major change in the 
mass range of Population III stars as halo masses increase to T^ir ~ 10^ K. This IS m 
mild disagreement with previous semianalytic results by Oh & Haiman (2002) and Mackey 
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et al. (2003). This is apparently because previous authors assumed that H2 coohng in 
"second generation" {Tyi^ > 10^ K) objects is unimportant, while we observe that a small 
H2 fraction, which forms in equilibrium with the UV background, can cool gas effectively to 
temperatures of ~ 1000 K. We do agree with these authors that the majority of primordial 
star formation will take place in objects with T^ir ~ 10"^ K, the so-called "Second generation" 
objects (Maclntyre, Santoro, & Thomas 2006) 

Our work appears to contradict the results of Greif & Bromm (2006), who use analytic 
and semianalytic techniques to follow the collapse of gas within halos of T^r > 10"^. They 
find that, in their model, gas cools to 8000 K via atomic line cooling, and then contracts 
nearly isothermally to high densities, allowing a molecular hydrogen fraction of Ih2 ~ 0.001 
to build up independent of local density and temperature. The gas in the center of the halo 
then cools and allows the gas to fragment on scales of ~ 20 M0. We do not see this mode 
of star formation, possibly because we never generate the significantly ionized halos with 
T> 10^ K upon which their scenario depends. Rather, we simulate halos that in the cases 
of J21 > 0.1 are built up to T^r ~ lO'^ K via mergers, and the gas is never ionized to a 
significant degree. One possible reason for the observed differences is that we ignore the 
formation of deuterium hydride and its effects on the cooling properties of primordial gas. It 
is possible that, in some contexts, the inclusion of HD may result in enhanced fragmentation. 
This will be examined in future studies. The scenario described by Oh & Haiman may still 
in fact be able to occur, and deserves further detailed numerical study to determine the fate 
of the collapsing gas and the implications of the (presumably) extremely high accretion rates 
onto the forming protostellar clouds. 

Our results are similar to those shown by Wise & Abel (2007), who perform similar sim- 
ulations examining the collapse of gas in cosmological halos in the presence of a photodissoci- 
ating background. They also find that, regardless of the strength of the UV background, the 
collapse of gas in the centers of the cosmological halos in question still occurs via H2 coohng. 
This holds even for quite extreme examples, such as when the original electron fraction is 
set to unphysically low levels (which should suppress H2 formation, at least temporarily). 

10. Summary 

In this paper we have performed a suite of high dynamical range (L/Ax ~ 5 x 10^) 
3D adaptive mesh cosmological simulations of the formation of Population III stars in a 
ACDM universe in the presence of a molecular hydrogen photodissocating ( "Lyman- Werner" ) 
ultraviolet background. The purpose of these calculations is to understand the effect that 
the soft ultraviolet background has on the evolution of the gas in collapsing cosmological 
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halos and to determine possible effects on the forming protostellar cloud. These simulations 
are all of a single cosmological realization, but with varied ultraviolet background strengths. 
Our principal results are as follows: 

1. Our calculations show that, as the flux of ultraviolet radiation in the Lyman- Werner 
band is increased. Population III star formation in a given cosmological halo is delayed 
to later times and, as a result, an increase in the virial mass of this halo at the onset of 
baryon collapse. This is in good agreement with previous work by Machacek et al. (2001) 
and Yoshida et al. (2003). 

2. We find that, contrary to previous work, the formation of primordial stars is never 
completely suppressed, regardless of the strength of the UV background. A small amount of 
molecular hydrogen always exists in cosmological halos, and allows cooling and halo collapse 
in gas which is bathed in a strong photodissociating background but has not reached lO'^ K. 
The previously-suggested mode of star formation in "second generation" halos, where the 
collapse of gas to high density is completely suppressed until cooling can take place via atomic 
hydrogen line transitions, is never observed in our calculations, which have photodissociating 
background strengths up to J21 = 1. 

3. Though the molecular hydrogen fraction in the centers of halos which are bathed in 
strong (J21 > 0.1) Lyman- Werner radiation is very small at early times {iH2 ~ 10~^ — 10^^), 
the gas at the center of the halo can still collapse due to the efficient cooling of molecular 
hydrogen at 2 — 5 x 10^ K, and to the extended time these halos require to collapse compared 
to halos in the presence of much smaller UV backgrounds. 

4. We observe that the estimated accretion rate onto the forming protostellar cloud 
varies strongly as a function of Jlw, with simulations that have a stronger ultraviolet back- 
ground having higher accretion rates. This is a function of the suppression of molecular 
hydrogen formation (and thus suppression of effective cooling) by the photodissociating 
background and of the higher virial temperatures of these halos when the gas finally col- 
lapses. 

5. Only a single clump is formed at the center of each collapsing cosmological halo, 

regardless of the strength of the photodissociating background, up to a baryon number 
density of n ~ 10^° cm~^. This implies that, as in the more commonly-studied mode of star 
formation, we will find only a single star per halo even in objects which are massive enough 
that T^ir > 10^ K. 
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